この文章で解説するのは
の4つです。ここまではおよそどのようなポイントクラウドデータでも 共通ですので、ここまでできるようになると良いでしょう。
import homcloud.interface as hc # HomCloudのインターフェス
import homcloud.plotly_3d as p3d # 3次元可視化用
import plotly.graph_objects as go # これも3次元可視化用
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
データを np.loadtxt で読みこみます。
pointcloud = np.loadtxt("pointcloud.txt")
まずは pointcloud の情報を調べてみましょう.
pointcloud.shape
np.min(pointcloud, axis=0), np.max(pointcloud, axis=0)
データは 1000x3 の配列です.これは3次元の点が1000個あることを意味します. またX,Y,Z座標の最大最小を見るとデータは$[0, 1]\times[0, 1]\times[0, 1]$という立方体に分布していることがわかります. この点を可視化してみます.
go.Figure(
data=[p3d.PointCloud(pointcloud, color="black")],
layout=dict(scene=dict(xaxis=dict(visible=False), yaxis=dict(visible=False), zaxis=dict(visible=False)))
)
可視化には plotly というライブラリを使っており,ブラウザ上のjupyter notebookで直接表示できます. マウスの左ボタンドラッグで回転,右ボタンドラッグで平行移動,ホイール回転で拡大縮小ができます.
dict(xaxis=dict(visible=False), yaxis=dict(visible=False), zaxis=dict(visible=False)) というのを略記するために p3d.SimpleScene という関数が使えます.
# 上と同じ
go.Figure(
data=[p3d.PointCloud(pointcloud, color="black")],
layout=dict(scene=p3d.SimpleScene())
)
データからパーシステント図を計算します。hc.PDList.from_alpha_filtrationという静的メソッドを使います。
save_boundary_map は後で optimal volume というものを計算するのに必要です。
hc.PDList.from_alpha_filtration(pointcloud,
save_to="pointcloud.pdgm",
save_boundary_map=True)
すると pointcloud.pdgm というファイルが生成されます。これが
パーシステント図の情報を収めたファイルです。
HomCloud 3.0 以降では .pdgm という拡張子を推奨しています.
pdlist = hc.PDList("pointcloud.pdgm")
このオブジェクトはすべての次数のパーシステント図を保持している(0次から2次まで)ので、1次のパーシステント図だけ取り出します。
すべての次数を保持しているので、このクラス自体の名前はPDListという名前です。
pd1 = pdlist.dth_diagram(1)
このメソッドの返り値は PD class のインスタンスです。
ヒストグラムを構築し、パーシステント図をプロットします。
histogramメソッドでヒストグラムを構築し、plotでそれをプロットします。
pd1.histogram().plot()
左下のほうに何かちょっと見えてあとはかなりぼんやりしています。というのは このデータは左下のほうにデータ(birth-death pair)が偏っているからです。 そこで色付けを log scale にしましょう。
pd1.histogram().plot(colorbar={"type": "log"})
基本的にはパーシステント図の点は対角線から離れるほど「意味のある」リング構造と対応し、 またY軸が大きい値になるほど大きなリング構造を表しています。 つまり(0.0025, 0.008)のあたりにある点がリング構造のなかで最もちゃんとしたリング構造を持っているものと対応していると言えそうです。
ここで注意しておくと、パーシステント図のX、Y軸は点に貼り付ける球の半径と対応していると
よくある教科書には書かれています。HomCloudでは実はこれは半径の2乗が使われています。
つまり$\sqrt{0.0025}=0.5$と$\sqrt{0.008} \approx 0.09$
が実際の半径になります。これは主には内部のアルゴリズムの都合に
よるものですが、各点に重みづけをしたときにはこのほうが自然に見えるという事情もあり、
HomCloudでは2乗の値が使われます。これを止めたいときはパーシステント図をhc.PDList.from_alpha_filtration
で計算するときに no_squared=True という引数を付けると
半径パラメータそのものがX、Y軸に現れます。
さて、次にbirth-death pairが集中している
左下の 0.0〜0.01 のあたりを拡大して調べてみましょう。デフォルトでは
ヒストグラムを構築するときにすべてのbirth-death pairが含まれるような
範囲を自動的に切り取ります。これを変更するには histogram メソッドのx_range引数(第1引数)を指定します。
pd1.histogram((0, 0.01)).plot(colorbar={"type": "log"})
ヒストグラムの解像度を変更するには histogram メソッドの
x_bins引数(第2引数)を指定します。デフォルトでは
128x128でヒストグラムを描きますが、もっと細かく256x256にしてみましょう。
pd1.histogram((0, 0.01), 256).plot(colorbar={"type": "log"})
これらの図は matplotlib の機能を使って描画されているため、
matplotlib.pyplot.savefigで保存することができます。以下のようにすると、pointcloud-pd1.pngに図が保存されます。
pd1.histogram((0, 0.01), 256).plot(colorbar={"type": "log"})
plt.savefig("pointcloud-pd1.png")
2次元のPDも同様に可視化しましょう.
pd2 = pdlist.dth_diagram(2)
pd2.histogram().plot(colorbar={"type": "log"})
こちらは death (Y軸) が 0.0125 より大きいあたりにぱらぱらと birth-death pair が分布しているのが特徴的です.
pd1 や pd2 が指すオブジェクトの births、deaths属性を調べることで、birth time、death timeがわかります。
pd1.births
pd1.deaths
そこで、この2つの差を見ることでlifetimeが計算できます。
pd1.deaths - pd1.births
lifetime のヒストグラムを表示してみましょう。
plt.hist(pd1.deaths - pd1.births, bins=100);
行末の ; としているのは plt.hist の返り値を無視するためのトリックです。大半の lifetime は非常に小さいことがわかります。実際対角線のそばに多くの birth-death pair がノイズ的な情報として分布していることはよくあります.
パーシステント図の個々の点は何らかの意味で元のポイントクラウドのリング構造、空隙構造と対応しているはずです。しかしそれがどのようなものであるのかを特定するのは実は そんなに簡単ではないです。このような解析を逆解析と呼びます。
ここでは HomCloud の強力な逆解析ツールである、optimal volumeを使いましょう。とりあえず (0.0025, 0.008) 付近にある birth-death pair の optimal volume を調べます。 optimal volume について詳しくは、大林の論文を参考にしてください。
pd.nearest_pair_to で (0.0025, 0.008) に一番近い birth-death pair を検索することができます。
pair = pd1.nearest_pair_to(0.0025, 0.008)
このオブジェクトは Pair クラスのオブジェクトで、birth timeやdeath timeといった情報を保持しています。
pair
以下のようにすると optimal volume が計算できます。
optimal_volume = pair.optimal_volume()
boundary_points メソッドで、対応するリング構造に含まれる点の座標が得られます。
optimal_volume.boundary_points()
このリング構造は39点からなっていることがわかります。boundaryメソッドを使うとどの点とどの点がつながっているか等もわかります。
このoptimal volumeを3次元可視化してみましょう。
go.Figure(data=[optimal_volume.to_plotly3d()], layout=dict(scene=p3d.SimpleScene()))
リングが表示されます。
次のようにするとポイントクラウドを重ねて表示することもできます.
go.Figure(data=[
optimal_volume.to_plotly3d(width=4, name="Optimal Volume"),
p3d.PointCloud(pointcloud, color="black", name="Pointcloud")
], layout=dict(scene=p3d.SimpleScene()))
to_plotly_meshを使うとループの内側の面を貼るような可視化も可能です.見やすい可視化を工夫してみてください.
fig.update_traces(opacity=0.5, selector=dict(name="Optimal Volume")) は面の部分の不透明度を50%にしています.
fig = go.Figure(data=[
optimal_volume.to_plotly3d(width=4, name="Optimal Volume outline"),
optimal_volume.to_plotly3d_mesh(color="red", name="Optimal Volume"),
p3d.PointCloud(pointcloud, color="black", name="Pointcloud")
], layout=dict(scene=p3d.SimpleScene()))
fig.update_traces(opacity=0.5, selector=dict(name="Optimal Volume"))
次のようにするとブラウザの別のタブで開きます.最後の行のfig.show(renderer="browser")が鍵です.
fig = go.Figure(data=[
optimal_volume.to_plotly3d(width=4, name="Optimal Volume outline"),
optimal_volume.to_plotly3d_mesh(color="red", name="Optimal Volume"),
p3d.PointCloud(pointcloud, color="black", name="Pointcloud")
], layout=dict(scene=p3d.SimpleScene()))
fig.update_traces(opacity=0.5, selector=dict(name="Optimal Volume"))
fig.show(renderer="browser")
name=...という引数で各オブジェクトに名前を付けたり,width=4という引数で線の太さを変更したりできます.
こんどは複数の点のoptimal volumeを同時に可視化してみます.birth < 0.004, death > 0.006 くらいの領域の birth-death pair を 以下のようにして取り出してきます.
pairs = [pair for pair in pd1.pairs() if pair.birth_time() < 0.004 and pair.death_time() > 0.006]
取り出してきた birth-death pair は以下のようにして数えることができます.
len(pairs)
では optimal volume を計算して,可視化してみましょう.54個もあるので計算には結構時間がかかります.そこで
cutoff_radiusというのを指定します.これは計算ターゲットが存在する領域の大きさを指定することで計算時間を
削減するための仕組みです.だいたいこの球の内側にあるだろう,という領域の半径を指定します.
この半径を予想するのは簡単ではありませんがここではデータ全体の座標が0〜1の間なので適当に1/3くらいの
0.3でいいかな,と予想してやりましょう.
%%time
optimal_volumes = [pair.optimal_volume(cutoff_radius=0.3) for pair in pairs]
go.Figure(data=[
v.to_plotly3d(width=4, name=f"({v.birth_time():.5f}, {v.death_time():.5f})") for v in optimal_volumes
] + [
p3d.PointCloud(pointcloud, color="black", name="Pointcloud")
], layout=dict(scene=p3d.SimpleScene()))
空洞のほうも同様に可視化してみましょう.birth > 0.0125 かつ death - birth > 0.0005 などがおもしろそうです.
pairs = [pair for pair in pd2.pairs() if pair.birth_time() > 0.0125 and pair.lifetime() > 0.0005]
optimal_volumes = [pair.optimal_volume() for pair in pairs]
go.Figure(data=[
v.to_plotly3d(width=4, name=f"({v.birth_time():.5f}, {v.death_time():.5f})") for v in optimal_volumes
] + [
p3d.PointCloud(pointcloud, color="black", name="Pointcloud")
], layout=dict(scene=p3d.SimpleScene()))
次のようにすると空隙部分の体積を可視化することもできます.update_tracesについては https://plotly.com/python/creating-and-updating-figures/ などを参考にしてください.
fig = go.Figure(data=[
v.to_plotly3d_mesh(name=f"({v.birth_time():.5f}, {v.death_time():.5f})") for v in optimal_volumes
] + [
p3d.PointCloud(pointcloud, color="black", name="Pointcloud")
], layout=dict(scene=p3d.SimpleScene()))
fig.update_traces(opacity=0.5, selector=dict(type="mesh3d"))